Bare Necessities

Author

Adam Kemberling

Published

November 11, 2024

Code 
# Libraries
library(tidyverse)
library(gmRi)
library(rnaturalearth)
library(scales)
library(sf)
library(gt)
library(patchwork)
library(lmerTest)
library(emmeans)
library(merTools)
library(tidyquant)
library(ggeffects)
library(performance)
library(gtsummary)
library(gt)
library(sizeSpectra)
library(ggdist)
library(pander)
library(ggh4x)
library(nlme)
library(ggimage)
conflicted::conflict_prefer("select", "dplyr")
conflicted::conflict_prefer("filter", "dplyr")
conflicted::conflicts_prefer(lmerTest::lmer)

# Processing functions
source(here::here("Code/R/Processing_Functions.R"))

# ggplot theme
theme_set(
  theme_gmri(
    axis.line.y = element_line(),
    axis.ticks.y = element_line(), 
    rect = element_rect(fill = "white", color = NA),
    panel.grid.major.y = element_blank(),
    strip.text.y = element_text(angle  = 0),
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    legend.position = "bottom"))



# labels for factor levels
area_levels <- c("GoM", "GB", "SNE", "MAB")
area_levels_long <- c("Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight")
area_levels_all <- c("Northeast Shelf", area_levels)
area_levels_long_all <- c("Northeast Shelf", area_levels_long)

# table to join for swapping shorthand for long-hand names
area_df <- data.frame(
  area = c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "Northeast Shelf"),
  survey_area = c("SS", "GoM", "GB", "SNE", "MAB", "Northeast Shelf"),
  area_titles = c("Scotian Shelf", "Gulf of Maine", "Georges Bank", "Southern New England", "Mid-Atlantic Bight", "Northeast Shelf"))


# Degree symbol
deg_c <- "\u00b0C"



# Function to get the Predictions
# Flag effect fits that are non-significant  ###
get_preds_and_trendsignif <- function(mod_x){
  modx_preds <- as.data.frame(
    # Model predictions
    ggpredict(
      mod_x, 
      terms = ~ est_year * area * season) ) %>% 
    mutate(
      area = factor(group, levels = area_levels_long_all),
      season = factor(facet, levels = c("Spring", "Fall")))
  
    # Just survey area and year
    modx_emtrend <- emtrends(
      object = mod_x,
      specs =  ~ area*season,
      var = "est_year") %>% 
      as_tibble() %>% 
      mutate(
        zero = 0,
        non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))
  
    # Preds with signif
    modx_preds %>% left_join(select(modx_emtrend, area, season, non_zero))
  
}

Storycrafting - Northeast Shelf Spectra

I think we need to step back a bit and find the main story again–which is about how size structure of the NES fish community has changed over time, whether that varies by region, and what covariates are associated with those changes.
The Aug draft you developed really focused on comparing regression models to consider drivers based on the fish community size spectrum as viewed from the perspective of the Wigley species set vs. the full species set.
However, I didn’t see figures of those spectra in that draft (maybe there are somewhere else and I’ve lost track). The recent slide deck really focuses on seasonal differences in the size spectra and influence of temperature. - KMills

Part 2: Relationships to External Factors

The aim of this section is to provide context and potentially attribution/correlations that may help explain trends/changes observed in the above section.

I will lead with a characterization of broad regional changes in temperature and fisheries landings totals.

Temperature and Size Structure

From what I see the size distribution along the shelf is responding consistently with temperature changes on shorter time scales than hypotheses that predict temperature related changes in growth.

The extent that the shape of the size distribution within each region follows some negative linear relationship to increasing temperature largely coincidental with the seasonal differences.

Seasonal movements contributing to this could maybe be explained by expectations from the MTE. But they also could be separate ecological processes like seasonal foraging or spawning opportunities. Either way, seasonal differences are large.

The consistency in the seasonal restructuring suggests to me that the size distribution is maybe less sensistive to long-term trends in temperature than other immediate factors. I say this noticing that the distribution is able to effectively track swings in temperature that happen within the year.

Code 
wigley_bmspectra_model_df %>% 
  #filter(survey_area == "Northeast Shelf") %>% 
  ggplot(aes(bot_temp, b)) +
  geom_point(aes(color = season)) +
  geom_smooth(method = "lm", color = "black") +
  facet_wrap(~survey_area, ncol = 1, labeller = label_wrap_gen(width = 8)) +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  theme(legend.position = "right") +
  labs(
    title = "Regional Temperature-Spectra Relationship\nReflects Seasonal Community Shift",
    subtitle = "Seasonal differences in spectra proportional to seasonal temperature range",
    shape = "Survey Area")

If we follow the temperature-spectra relationship within a particular season e.g. GOM-Spring, the relationship is much weaker and not-significant in most cases.

I think there is an important question about temperature’s role, but I don’t think that local warming is having much impact on the community size distribution. If it is, it is secondary to other size distribution shifts happening at shorter time scales and across spatial scales.

This raises questions about how seasonal size distribution changes are being achieved.

From the decadal variability work we have reason to believe that individual species are both not tracking temperature changes in space at the rate they are ocurring AND that changes in size at age are happening. These are signals at the species level indicating that temperature effects on growth are ocurring.

There would need to be some compensatory mechanisms present to offset this across the broader community.

Framing “Warming” Using BT Anomalies

I want to avoid conflating that a relationship between spectra and bottom temperature means that “warming” has altered the size spectra. There are known relationships between inter-specific and intra-specific body size across temperature ranges and across latitudinal ranges.

  • Bergmann’s Rule: Warmer climates are often inhabited by smaller-bodied species
  • James Rule: At an intraspecific level, populations living in warmer conditions often reach smaller maximum body sizes

The metabolic theory of ecology also supports/provides a basis for the above rules, as well as seasonal redistribution related to body size and temperature.

  • Metabolic Theory of Ecology: Metabolic theory predicts how metabolic rate, by setting the rates of resource uptake from the environment and resource allocation to survival, growth, and reproduction, controls ecological processes at all levels of organization from individuals to the biosphere.

Metabolic Efficiency and Seasonal Movement: MTE suggests that metabolic rates, which are influenced by body size and temperature, are central to the ecological behavior of organisms.

Larger individuals, which naturally have lower per-unit mass metabolic rates, may move to cooler areas seasonally to optimize their metabolic function and reduce the energetic costs of thermoregulation in warmer periods.

Code 
wigley_bmspectra_model_df %>% 
  filter(survey_area != "Northeast Shelf") %>% 
  ggplot(aes(bot_temp, b, color = season)) +
  geom_point(aes(shape = survey_area)) +
  geom_smooth(method = "lm", color = "black") +
  ggforce::geom_mark_ellipse(aes(fill = season), alpha = 0.1) +
  scale_color_gmri() +
  scale_fill_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  guides(
    shape = guide_legend(nrow = 2),
    fill = guide_legend(nrow = 2),
    color = guide_legend(nrow = 2)) +
  theme(legend.box = "horizontal") +
  labs(
    title = "Bergmann's Rule + MTE",
    subtitle = "Broad relationship between temperature and size distribution,\nbolstered by seasonal temperature-consistent movements",
    shape = "Survey Area",
    fill = "Season",
    color = "Season",
    x = "Bottom Temperature")

The observation that body size distribution broadly follows some temperature relationship is a confirmation of these broadly observed rules, but I would argue it does not convey that warming has shifted a local climate in some direction.

Since we know that individuals migrate and may use the full water column or occupy different habitats seasonally, we need to treat seasonal communities separately to avoid an apples to oranges comparison of the community size distribution.

Using these same season + region blocks, we can characterize the typical bottom temperature for these location and seasonal combinations using the long-term average, and quantify warming using departures from this long-term average.

This also resolves a separate issue when modeling with absolute temperatures, where temperature and season*area units are colinear.

Local Bottom Temperature Anomalies Model

This is the model I would suggest gets the most at “Warming” impacting regional size distributions:

Fixed effects:

  • Season
  • Survey area
  • Bottom temperature anomaly

Error structure: AR1 autocorrelative errors

Code 
# Just drop the random effect part for now
anom_model_ar <- nlme::gls(
  b ~ area * season * scale(bot_temp_anom), 
  data = wigley_bmspectra_model_df, 
  correlation = corAR1(form = ~ est_year | area/season)
)


# # confidence intervals on phi
# intervals(anom_model_ar)$corStruct


# # check the model
# check_model(anom_model_ar)
# plot(check_collinearity(anom_model_ar))
# plot(
#   check_collinearity(
#     nlme::gls(
#       b ~ area + season + scale(bot_temp_anom), #+ scale(log(total_weight_lb)), 
#       data = wigley_bmspectra_model_df, 
#       correlation = corAR1(form = ~ est_year | area/season)
# ))
#   )


# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    anom_model_ar, 
    terms = ~ bot_temp_anom * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####

# Flag effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = anom_model_ar,
  specs =  ~ area * season,
  var = "bot_temp_anom") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  group_by(season, area) %>%
  summarise(
    min_temp = min(bot_temp_anom)-2,
    max_temp = max(bot_temp_anom)+2,
    .groups = "drop")



# Plot over observed data
local_btemp_anoms <- mod2_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_temp) == F,
         (x > max_temp) == F) %>%
  left_join(select(mod2_emtrend, area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wigley_bmspectra_model_df,
    aes(bot_temp_anom, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(y = "Mass distribution Exponent",
       x = "Local Seasonal Bottom Temperature Anomaly")
local_btemp_anoms

Based on this model I would say there is evidence of a warming effect in: GOM-Spring, GB-Spring, GB-Fall, & SNE-Fall.

Temperature Model

We were originally using absolute temperatures, which are scaled, and I was curious with the AR error structure how different this would be.

I just had to see if the same model is better performing or leads to different results when it knows absolute temperatures.

Code 
# Just drop the random effect part for now
temp_model_ar <- nlme::gls(
  b ~ area * season * scale(bot_temp), 
  data = wigley_bmspectra_model_df, 
  correlation = corAR1(form = ~ est_year | area/season)
)


# # confidence intervals on phi
# intervals(temp_model_ar)$corStruct
# 
# # check the model
# check_model(temp_model_ar)

# # Collinearity without interactions
# plot(
#   check_collinearity(
#      nlme::gls(
#         b ~ area + season + scale(bot_temp), 
#         data = wigley_bmspectra_model_df, 
#         correlation = corAR1(form = ~ est_year | area/season))
# ))


# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    temp_model_ar, 
    terms = ~ bot_temp * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####

# Flag effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = temp_model_ar,
  specs =  ~ area * season,
  var = "bot_temp") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  group_by(season, area) %>%
  summarise(
    min_temp = min(bot_temp)-2,
    max_temp = max(bot_temp)+2,
    .groups = "drop")



# Plot over observed data
local_btemp_results <- mod2_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_temp) == F,
         (x > max_temp) == F) %>%
  left_join(select(mod2_emtrend, area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wigley_bmspectra_model_df,
    aes(bot_temp, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(y = "Mass distribution Exponent",
       x = "Seasonal Bottom Temperature")
local_btemp_results

Verdict: Temperature itself performs better, and is simpler to produce/explain.

Code 
data.frame(
  "Model" = c("Temperature Model", "Anomalies Model"),
  "AIC" = c(AIC(temp_model_ar), AIC(anom_model_ar))) %>% 
  gt()
Model AIC
Temperature Model -1206.981
Anomalies Model -1185.086

Total Landings

The landings information is the newer GARFO data obtained from Andy Beet. It was delivered as containing only finfish, I have removed freshwater species and blue crabs. It should largely reflect marine finfish species at this point.

Code 
# From Andy Beet - Crabs removed
landings_annual <- read_csv(here::here("Data/processed/BEET_GARFO_regional_finfish_landings.csv")) %>% 
  select(survey_area,
         est_year = year,
         total_weight_lb = total_live_lb)

shelf_landings_annual <-  mutate(landings_annual, survey_area = "Northeast Shelf") %>% 
    group_by(est_year, survey_area) %>% 
    summarise(total_weight_lb = sum(total_weight_lb, na.rm = T))



# Plot annual totals
bind_rows(list(
  shelf_landings_annual,
  landings_annual)) %>% 
  mutate(survey_area = factor(survey_area, levels = area_levels_long_all)) %>% 
  filter(est_year <= 2019) %>% 
  ggplot(aes(est_year, total_weight_lb/1e6)) +
  geom_col(alpha = 0.8) +
  geom_ma(
    aes(color = "5-Year Moving Average"), 
    size = 1, n = 5, ma_fun = SMA, key_glyph = "timeseries", linetype = 1) +
  scale_color_manual(values = "orange") +
  facet_wrap(~survey_area, ncol = 1, scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_y_continuous(labels = label_number(suffix = "M")) +
  labs(y = "Total Landings (live lb.)", x = "Year", color = "Moving Average")

Landings Autoregressive model

I wanted to follow on the bottom temperature anomalies model here for landings:

Code 
# Just drop the random effect part for now
f_model_ar <- nlme::gls(
  b ~ area * season + area * scale(log(total_weight_lb)), 
  data = drop_na(wigley_bmspectra_model_df, total_weight_lb) %>% 
    dplyr::filter(est_year > 1978),
  correlation = corAR1(form = ~ est_year | area/season)
)



# Clean up the plot:
# Plot the predictions over data
mod2_preds <- as.data.frame(
  ggpredict(
    f_model_ar, 
    terms = ~ total_weight_lb * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####
trend_df <- emtrends(
  object = f_model_ar,
  ~area,
  var = "total_weight_lb",
  mode = "appx-satterthwaite")


# Flag effect fits that are non-significant  ###
mod2_emtrend <- emtrends(
  object = f_model_ar,
  specs =  ~ area,
  var = "total_weight_lb",
  mode = "appx-satterthwaite") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  group_by(season, area) %>%
  summarise(
    min_f = min(total_weight_lb)-2,
    max_f = max(total_weight_lb)+2,
    .groups = "drop")



# Plot over observed data
landings_ar_plot <- mod2_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_f) == F,
         (x > max_f) == F) %>%
  left_join(select(mod2_emtrend, area,non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = wigley_bmspectra_model_df,
    aes(total_weight_lb, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free") +
  scale_color_gmri() +
  scale_x_log10(labels = label_log(base = 10), limits = 10^c(0,10)) +
  labs(y = "Exponent of ISD",
       title = "Total Landings (lb.)",
       x = "Local Seasonal Bottom Temperature Anomaly")
landings_ar_plot

Code 
# temp model is better
# AIC(temp_model_ar)
# AIC(f_model_ar)

This is our old friend “higher landings when spectra was shallower” that we’ve seen before.

Anomalies and Landings Together

If we’re using local/seasonal temperature anomalies we have less of an issue with singular fits for the region*season interaction and can more lucidly model with landings simultaneously.

Code 
# Just drop the random effect part for now
full_model_ar <- nlme::gls(
  b ~ area * season * scale(bot_temp_anom) + area * scale(log(total_weight_lb)), 
  data = drop_na(wigley_bmspectra_model_df, total_weight_lb) %>% 
    dplyr::filter(area != "Northeast Shelf"),
  correlation = corAR1(form = ~ est_year | area/season))

# # Collinearity without interactions
# plot(
#   check_collinearity(
#     nlme::gls(
#         b ~ area + season + scale(bot_temp_anom) +  scale(log(total_weight_lb)),
#         data = drop_na(wigley_bmspectra_model_df, total_weight_lb) %>% 
#                  filter(area != "Northeast Shelf"),
#         correlation = corAR1(form = ~ est_year | survey_area/season)
#       )
#     )
#   )


# # confidence intervals on phi
# intervals(full_model_ar)$corStruct


# temp model is better
# AIC(temp_model_ar)
# AIC(full_model_ar)
Code 
# Clean up the plot:
# Plot the predictions over data
temp_preds <- as.data.frame(
  ggpredict(
    full_model_ar, 
    terms = ~ bot_temp_anom * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####

# Flag effect fits that are non-significant  ###
temp_emtrend <- emtrends(
  object = full_model_ar,
  specs =  ~ area * season,
  mode = "appx-satterthwaite",
  var = "bot_temp_anom") %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  filter(area != "Northeast Shelf") %>% 
  group_by(season, area) %>%
  summarise(
    min_temp = min(bot_temp_anom)-2,
    max_temp = max(bot_temp_anom)+2,
    .groups = "drop")



# Plot over observed data
local_btemp_anoms <- temp_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_temp) == F,
         (x > max_temp) == F) %>%
  left_join(select(temp_emtrend, area, season, non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = filter(wigley_bmspectra_model_df, area!= "Northeast Shelf"),
    aes(bot_temp_anom, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_color_gmri() +
  scale_x_continuous(labels = label_number(suffix = deg_c)) +
  labs(y = "Mass distribution Exponent",
       x = "Local Seasonal Bottom Temperature Anomaly")
local_btemp_anoms

Code 
# Clean up the plot:
# Plot the predictions over data
f_preds <- as.data.frame(
  ggpredict(
    full_model_ar, 
    terms = ~ total_weight_lb * area * season))   %>% 
  mutate(
    area = factor(group, levels = area_levels_long_all),
    season = factor(facet, levels = c("Spring", "Fall")))



#### Trend Posthoc  ####

# Flag effect fits that are non-significant  ###
f_emtrend <- emtrends(
  object = full_model_ar,
  specs =  ~ area,
  var = "total_weight_lb",
  mode = "appx-satterthwaite"
  ) %>%
  as_tibble() %>%
  mutate(
    zero = 0,
    non_zero = if_else(between(zero, lower.CL, upper.CL), F, T))


# Limit temperature plotting range
actual_range <- wigley_bmspectra_model_df %>%
  group_by(season, area) %>%
  summarise(
    min_f = min(total_weight_lb)-2,
    max_f = max(total_weight_lb)+2,
    .groups = "drop")



# Plot over observed data
landings_ar_plot <- f_preds %>% 
  left_join(actual_range) %>%
  filter((x < min_f) == F,
         (x > max_f) == F) %>%
  left_join(select(f_emtrend, area,non_zero)) %>%
  filter(non_zero) %>%
  ggplot() +
  geom_ribbon(aes(x, ymin = conf.low, ymax = conf.high, group = season), alpha = 0.1) +
  geom_line(
    aes(x, predicted, color = season), 
    linewidth = 1) +
  geom_point(
    data = filter(wigley_bmspectra_model_df, area != "Northeast Shelf"),
    aes(total_weight_lb, b, color = season),
    alpha = 0.4,
    size = 1) +
  facet_grid(area~., scales = "free", labeller = label_wrap_gen(width = 8)) +
  scale_color_gmri() +
  scale_x_log10(labels = label_log(base = 10), limits = 10^c(0,10)) +
  labs(y = "Exponent of ISD",
       title = "Total Landings (lb.)",
       x = "Local Seasonal Bottom Temperature Anomaly")
landings_ar_plot

Zooplankton Community Indices

The ecodata package has changed the zooplankton indices around that it carries. It now contains “regime” metrics for 26 of the main zooplankton taxa.

https://noaa-edab.github.io/tech-doc/zoo_abundance_anom.html

Code 
zoo_regimes <- ecodata::zoo_regime
# zoo_anoms <- ecodata::zoo_abundance_anom
zoo_regimes %>% 
  filter(str_detect(Var, "calfin|chaeto|pseudo|oith")) %>% 
  mutate(EPU = factor(EPU, levels = c("GOM", "GB", "MAB"))) %>% 
ggplot(aes(Time, Value)) +
  geom_point(
    size = 0.8, alpha = 0.5) +
  geom_ma(
    aes(linetype = "5-Year Moving Average"),
    size = 1, n = 5, ma_fun = SMA, key_glyph = "timeseries") +
  #scale_color_gmri() +
  facet_grid(EPU~Var) +
  labs()

Driver Correlations